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Abstract. We study the flow of a shear-thinning, chemically- reacting fluid that could be used 
to model the flow of the synovial fluid. The actual geometry where the flow of the synovial fluid 
takes place is very complicated, and therefore the governing equations are not amenable to simple 
• mathematical analysis. In order to understand the response of the model, we choose to study 

the flow in a simple geometry. While the flow domain is not a geometry relevant to the flow 
of the synovial fluid in the human body it yet provides a flow which can be used to assess the 
\ efficacy of different models that have been proposed to describe synovial fluids. We study the 

flow in the annular region between two cylinders, one of which is undergoing unsteady oscillations 
about their common axis, in order to understand the quintessential behavioral characteristics of 
the synovial fluid. We use the three models suggested by Hron et al. [ J. Hron, J. Malek, P. 
Pustejovska, K. R. Rajagopal, On concentration dependent shear-thinning behavior in modeling 
of synovial fluid flow, Adv. in Tribol.(In Press).] to study the problem, by appealing to a 
semi-inverse method. The assumed structure for the velocity field automatically satisfies the 
constraint of incompressibility, and the balance of linear momentum is solved together with a 
| convection-diffusion equation. The results are compared to those associated with the Newtonian 

model. We also study the case in which an external pressure gradient is applied along the axis 
of the cylindrical annulus. 

> 

in 

r — . 

. 1. Introduction 

There are several biological fluids as well as chemically-reacting polymeric liquids wherein 
q the response characteristics of the fluids change due to the extent of the chemical reaction. In 
i ■ a biological material, such as blood, which in a homogenized sense is modeled as a fluid, the 
, number of constituents are manifold: plasma, red and white blood cells, platelets, and a plethora 
of proteins; moreover, there are literally dozens of biochemical reactions that take place in order 
to maintain the blood in a delicate state of balance. For such complex fluids, the modeling 
advocated in this paper would be inappropriate. However, there are other biological liquids such 
as the synovial fluid which are less complex in structure in that there are fewer constituents and 
there are far fewer biochemical reactions. The present study addresses such biological fluids or 
polymeric liquids in which few or preferably one dominant chemical reaction takes place. 

Synovial fluid is a mixture that responds in a viscoelastic fashion (see [I], [2], [3], and [1]). 
That the properties of the fluid change due to pathological conditions has been established by 
Gomez and Thurston [5]. When the synovial fluid is flowing under conditions where there are 
no instantaneous inputs, it essentially behaves as a generalized viscous fluid, namely a Stokesian 
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fluid, and it is only when it is subject to instantaneous inputs does its viscoelastic nature manifest 
itself. Thus, the study of the flow of a generalized viscous fluid has some bearing on the flows of 
the synovial fluid. Recently, Hron et al. [6] studied the flow of three different fluid models which 
could be used as a model for the synovial fluid. All the models that they used fall into the category 
of generalized viscous fluids. One of them is a shear-thinning fluid model in which the power-law 
exponent is dependent on the concentration which leads to very interesting issues concerning its 
mathematical features. It is important to assess which of these three models studied by Hron et 
al. [6] is better suited to model the synovial fluid. We seek a better understanding of the models; 
especially with regard to determining if there are significant differences in the results between 
the models when the flow takes place in a totally different geometry. We consider an idealized 
geometry which allows us to use a semi-inverse method to study the problem while Hron et al. 
[6] use a finite element approach. However, unlike the geometry which Hron et al. |6] used, 
namely the flow in a rectangular region, here we consider a geometry wherein the boundary has 
a curvature which presents a new feature with regard to the flow domain in which to assess the 
characteristics of the model. 

Synovial fluid is a mixture of plasma that is free of large proteins but containing polysaccharide 
molecules called Hyaluronan (HA). A detailed description of the composition and structure of 
the synovial fluid can be found in [7]. The presence of HA makes the fluid non- Newtonian, and 
the response of the synovial fluid changes with the concentration of HA. However, in view of 
the mass/volume fraction of the HA in the synovial fluid being small, we can approximate it 
essentially as a coexisting mixture of the plasma (which will henceforth be identified as the fluid) 
and a reactant (which is essentially HA). Let us consider a fluid whose property changes due to its 
reacting with a chemical agent whose mass fraction is significantly smaller than that of the fluid 
and furthermore the diffusion velocity of the chemical agent with respect to the fluid is far less 
than that of the fluid velocity. This implies that the relative velocity between the fluid and the 
diffusant is not significant, and the diffusant is essentially carried by the fluid. It is important to 
recognize that while the diffusion velocity with respect to the fluid is small, it cannot be ignored, 
for otherwise we would not have the diffusion taking place within the fluid. Essentially what we 
have, which is the counterpart to the diffusion of a fluid through a porous solid, is the diffusion of 
one fluid through another fluid that is flowing with a significantly higher velocity with respect to 
the relative velocity between the fluid and the diffusant. With respect to the velocity of the fluid, 
the diffusion velocity of the diffusant can be ignored. However, we do have a term that takes into 
account the changes in the concentration of the diffusant due to concentration gradients. The 
properties of the fluid change due to the concentration of the reactant, as the fluid chemically 
reacts with the diffusant. The study that is carried out is thus relevant not just to the flow 
of the synovial fluid, but also to polymeric fluids that are chemically-reacting with a diffusing 
substance, in our case (HA). Thus the approach used by Bridges et al. [8] can be used to study 
the flow in question, which is the methodology adopted by Hron et al. [6]. 

We study the flow of this reacting fluid mixture in the annular region between two coaxial 
cylinders, one of which is rotating; there is also flow due to the presence of a pressure gradient. 
The organization of the paper is as follows. The preliminaries and kinematics required for 
this paper are discussed first in section (j2J). The models chosen to characterize the response 
of the synovial fluid, and the corresponding governing equations are described in section ([3]) 
and section (j3J), respectively. In the following section, we describe the initial boundary value 
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problem, the non-dimensionalization scheme, and the solution procedure in detail. The results 
are discussed in the final section. 

2. Kinematics 

Let B denote an abstract body, and let kr and n t represent the mappings, referred to as placers, 
which take the abstract body into the reference configuration Kr{B) and current configuration 
n t (B), respectively. We shall choose to use the initial configuration as the reference configuration 
for the sake of convenience. Furthermore, we shall assume that there is a one-parameter family 
of placers which can be associated with the abstract motion of the body. We shall associate 
the parameter with time t, t e Wq. In view of this assumption, we can define a mapping 
X KR ■ kr{B) xRq £ such that: 

x = x K/{ (X,f) where X e k r (B) and x e K t (B). (1) 

The velocity of the particle X is defined through the relation: 

v<(Jl (X ) = % a (X,t):=^ J (2) 

where —— represents the material time derivative, and the deformation gradient is defined as 
at 

follows: 

F Kfl (X,t) = A[ Xftjj (X,t)] : = Grad[x Kfl (X,t)]. (3) 

We shall assume that det(F Kfl ) > 0. Furthermore, we assume that the motion is sufficiently 
smooth so that all of the derivative operations that follow have meaning. Since the motion is 
presumed to be one-to-one for each instant of time, the mapping (1) can be inverted so that: 

X = (4) 

and therefore any physical quantity of interest q can be expressed in the various forms: 

<; = q(X, t) = <T(x^(x, *),*) = <f(x, t). (5) 

We shall primarily be interested in stating the functions in terms of x, i.e., we wish to study the 
mechanics of the reacting fluid mixture from an Eulerian perspective. Furthermore, we shall use 
lower case letters to represent differential operators such as div(-), grad(-), etc, which are carried 
out with respect to x, and this is consistent with using capital letters to signify that the derivative 
is taken with respect to X, as done previously with the definition of the deformation gradient 
(Lagrangean gradient of the motion). In addition, we shall henceforth omit the subscript Kr. 
It follows that 

|[det(F)] = det(F)div(v), (6) 

and 

L := grad(v) = FF~\ (7) 

where ( ) denotes the material time derivative. Now, if the body is incompressible and is 
therefore only capable of undergoing isochoric motions, then det(F) = 1 for each instant of time, 
and the relationship (6) reduces to 

div(v) = 0. (8) 
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For the purpose of this work, we shall find it convenient to define the Rivlin-Ericksen tensors 
(see [9]). The first Rivlin-Ericksen tensor has the following representation: 

Ax = L + L T = 2D, (9) 

where D is the symmetric part of the velocity gradient, i.e., the stretching tensor. The higher 
order Rivlin-Ericksen tensors are defined through the recurrence relationship: 

A* = ^=i + A fc _ 1 L + L T A fc _ 1 , k>2. (10) 

We shall employ a constrained mixture theory to study the behavior of synovial fluid which 
shear-thins or thickens due to the chemical reactions which take place. While the synovial fluid is 
comprised of numerous constituents: plasma, HA, fats, and leukocytes, the rheological properties 
of the mixture primarily depend upon the amount of plasma and HA. Thus, it is assumed that at 
each point within the mixture, the plasma (fluid) coexists with the HA (reactant), and the flow 
is constrained so these constituents move together. Henceforth, when referring to the "fluid" we 
shall mean the plasma, and we use the term "reactant" to stand for HA. We shall account for 
thinning and thickening by allowing the viscosity of the fluid to depend upon the shear-rate and 
the concentration, which changes by virtue of the chemical reactions that take place within the 
mixture. Here we define the concentration in the following manner: 

c=^^, (11) 

Qr + Qf 

where the quantities Qf and g r represent the densities of the fluid and the coexisting reactant, 
respectively. We shall assume that the fluid is incompressible, and therefore only isochoric 
motions are possible. 



3. Constitutive Assumptions 

As mentioned in the introduction, the synovial fluid is a chemically-reacting, non-linearly vis- 
coelastic, incompressible fluid mixture which possesses instantaneous elasticity. In fact, numerous 
experiments have been carried out which illustrate these facts (pQ, [2], [I], [5], and j3]). However, 
if the fluid is not subject to instantaneous inputs, and if the inputs are sufficiently slow, then the 
fluid can be approximated as an isotropic and incompressible fluid of the differential type. The 
fluids of the differential type are incapable of instantaneous elastic response. The Cauchy stress 
for such a fluid has the following representation (see Truesdell and Noll [ID]): 

T = -7rl + g(A lj A 2 ,...,A fc ), (12) 

where — nl denotes the spherical stress due to the incompressibility constraint. Fluids that can 
be represented by this model are usually referred to as incompressible Rivlin-Ericksen fluids of 
complexity k, and these fluid models are frame-indifferent as the tensors A^ are frame- indifferent. 
For a more general class of fluids of the differential type of class (m, k), see Dunn and Rajagopal 
[TT] . From the standard representation theorems, for an isotropic, tensor- valued function it 
follows that the Cauchy stress for an incompressible Rivlin-Ericksen fluid of complexity 1 has 
the following form: 

T = -tt1 + 1 A 1 + 2 A2, (13) 

where <pi and fa are each functions of the invariants tr(A^) and tr(Af). We are interested in a 
very special sub-class of (13) that is usually referred to as power-law fluids and such fluids can be 
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used to describe shear-thinning a characteristic of the synovial fluid. For such fluids, the Cauchy 
stress reduces to: 

T = -ttI + At (A x ) Ax, (14a) 
and we shall assume that the apparent viscosity jji a has the following structure: 



£i a (Al) = jJiQ 



1 + 7tr(A?) 



(14b) 



where hq, 7, n are each constant. These fluids are capable of characterizing shear-thinning and 
shear-thickening phenomena, but cannot describe normal stress differences in simple shear flows. 
Furthermore, the model does not exhibit stress-relaxation effects nor does it possess the capability 
to describe instantaneous elasticity; the model is capable of characterizing certain types of non- 
linear creep, however. If the shear-index n < 0, then this model characterizes the behavior of a 
shear-thinning fluid while if n > 0, it characterizes the behavior of a shear-thickening fluid. If 
n — 0, then this model reduces to the incompressible, linearly viscous fluid model. The quantity 
/i is called the zero shear-rate viscosity and is defined as: 

Ho = lim (15) 

where k is the shear-rate in a simple shear flow. The Lagrange multiplier n, in general, is a 
function of x and t, and since tr(Ai) = 2 div(v) = 0, n is equal to the mean normal stress. This 
is not the case with most incompressible, non-linear fluids as the trace of the extra stress is, in 
general, non-zero. 

We shall assume that the fluid is initially of uniform temperature and is constrained to undergo 
only isothermal processes^. Consequently, the rate of dissipation E reduces to: 

S = T D, (16) 

where the scalar product between the Cauchy stress and the stretching tensor is often referred 
to as the stress-power. The 2 nd law of thermodynamics requires that the rate of dissipation, and 
in this case the stress-power, to be non-negative. If the stress-power is to be non-negative for 
all n, then we can conclude that 7^0. These fluids have been used in modeling the flows of 
polymeric liquids [12], and the mathematical properties pertaining to existence and uniqueness 
of the equations governing the flows of a fluid characterized by (14) have been analyzed in great 
detail (see [13]), and the stability of such flows has been explored by Malek et al. [H]. Finally, 
we note that such fluid models can also be developed within the thermodynamic framework 
recently proposed by Rajagopal and Srinivasa [15] . In fact, if it turns out that the temperature or 
viscoelasticity of the fluid should be taken into consideration, then the thermodynamic framework 
used in [13] is much more appropriate and leads to a much richer class of models. 

As stated previously, we are interested in capturing the response of the synovial fluid mixture 
which not only shear-thins, but also chemically-thickens, and therefore we are interested in 
generalizations of the model (14). The model is generalized to account for chemical-thickening 
by allowing the apparent viscosity to also depend upon the concentration. There are a number 
of ways one could go about doing this; we shall be interested in three different models. The first 
model is due to Lai et al. [16] and Laurent et al. [17] . and it allows the zero shear-rate viscosity 



"The human body temperature continually fluctuates, and both the magnitude and manner in which these fluc- 
tuations take place can depend on a number of variables, both internal and external. 
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to change with concentration. In particular, the zero shear-rate viscosity takes the following 
form: 

/i (c) := /I e" c , (17) 

where ~p is the viscosity of the plasma, a is a constant. The other models we shall consider 
allow for a much more intricate connection between the stretching and concentration, and these 
models are due to Hron et al. |6|. For these models, which we shall designate as model 2a and 
2b, respectively, we allow the power-law exponent n that appears in f!14bj) to depend on the 
concentration in the specific way summarized as follows: 



nic 



2 v 



1) 



and 



n(c) :- 



a 



1 



ac 2 + 1 



(19) 



respectively, where a and a are each constant. Also, for these models wherein the shear index 
depends upon the concentration, the zero shear-rate is constant, i.e., fiQ = JI . Hron et al. [6] 
show that the structure of these models allows for a better fit to relative viscosity (=-) / shear- 
rate data, for a variety of concentrations than the models discussed in [16] and [17]. Table 1 
summarizes the constants used in each of the models. Thus, the general form of the models that 
we shall consider has the representation: 



-7Tl + //(c,Ai)Ai, 



where 



fi(c, Ai 



£ + 7tr(A?; 



n(c) 



(20a) 



(20b) 



Note that the constant (3 has also been incorporated into the model (20) which allows for a better 
fit to the test data [6]. A sufficient condition for the rate of dissipation to be non-negative, for all 
n(c), is to require both 7^0 and (5 ^ 0, and the parameter values obtained from the regression 
analyses obviously agree with this constraint. 



Table 1. Values for the Parameters used in each Model 



Parameter 


Model 1 


Model 2a 


Model 2b 


a 


21.3 


3.3 


31.0 


P 


1 


7.1 x lO" 9 


1.3 x 10~ 8 


7 t 


6.96 
4 


5.8xl(T 8 
4 


8.5xl0" 8 


a 






0.44 


n 


-0.28 







Note that 7 defined in model (20) is j times that 
used by Lai et al. [TB] , Laurent et al. [T7] , and Hron 
et al. [6] as their models are written in terms of D 
and the models herein are given in terms of Ai . 
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4. Governing Equations 

As stated previously, we employ the framework of constrained mixtures to model the synovial 
fluid mixture, and the constituents are therefore constrained to move together. Thus, we need 
not concern ourselves with the equations of motion for the reactant, just the equations which 
govern the flow of the fluid and the reactions which take place therein. We shall assume that 
the plasma is incompressible, and therefore we must ensure that the constraint (8) is satisfied. 
Thus, the balance of mass for the fluid: 

^ + ^div(v) = 0, (21) 

reduces to — ^ = 0. This implies that the fluid density is constant at each material point X. 
at 

However, since the fluid (plasma) is assumed to be homogeneous, it follows that Qf is the same 
constant everywhere. On substituting the constitutive relation (20) into the balance of linear 
momentum for the fluid: 

<9v 



Qf 

we obtain: 



+ grad(v)v 



div(T) + Q f b, (22) 



Qf 



<9v 

- + div( Vi 



-grad(7r) + div[/i(c, Ai)Ai]. (23) 



Note that we have utilized the condition (8) and neglected the body force Qfh in deriving equation 
(23). We shall assume that the chemical reactions which take place within the mixture are 
governed by the convection-diffusion equation: 

dc 

— + div(cv) = -div(p), (24) 

where p is a flux vector that is related to the reactions that take place within the synovial fluid. 
We shall assume that the flux vector p is given by a constitutive relation that is similar to that 
used in Fick's assumption, i.e., 

p = -Kgrad(c), (25) 

where K is, in general, a tensor-valued function of the concentration and the first Rivlin-Ericksen 
tensor. Therefore, the convection-diffusion equation reduces to: 

dc 

— + grad(c) • v = div [K(c, Ai)grad(c)] (26a) 

after enforcing the constraint (8). In fact, the diffusivity could depend on A*, k ^ 2, but we 
shall not consider such a possibility here. Experimental data suggests that, for the synovial fluid, 
the dependence on the concentration and the first Rivlin-Ericksen tensor is mild, and therefore 
we shall assume the diffusivity has the following form: 

K(c,Ai)=D c l, (26b) 

where D c is a constant. 

Thus, in general, the isochoricity constraint (8), the balance of linear momentum (23), and 
the convection-diffusion equation (26) are coupled, and one must solve them together for the 
velocity v, the pressure n, and the concentration c on a complicated, time-varying domain. Hron 
et al. [S] used a finite element approach to study the problem, and their approach can be used 
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in complicated flow geometries. In this study, our interest lies in carrying out the analysis in a 
different but simpler geometry in view of the difficulties involved with solving the problem in a 
complex geometry, especially one in which the curvature effects come into play. Thus, we shall 
consider a relatively simple initial boundary value problem. 

5. An Initial Boundary Value Problem 

5.1. Assumptions for the Velocity Field. As with any semi-inverse approach, one makes an 
attempt to reduce the complexity of the governing equations so that they become more tractable 
in a coordinate system appropriate for the problem of interest. The aim is to glean some insight 
into the full set of equations by studying the solution of the resulting simplified equations, usually 
on an infinite domain. For the problem under consideration, one such approach would be to study 
the plane flow wherein the velocity field has the form: 

v = v{u,^,t)e v + w{u,^,i)e^, (27) 

in a confocal-ellipsoidal coordinate system Of course, in such a coordinate system the 

lines of constant C form ellipsoids, while the lines of constant v and £ form hyperboloid of one 
sheet and two sheets, respectively; these surfaces could characterize the articular membrane 
which surrounds the bone and the synovial membrane of the joint. The drawback of such an 
approach is that as the appendage moves, the bone rotates and the coordinate system is no 
longer useful as the boundaries do not coincide with lines of constant hyperboloids. 

Another approach is to study the response of the reacting fluid mixture that fills to the gap 
between two spheres and allow the outer sphere to rotate while the inner sphere remains fixed. 
One could utilize the similarity transformation suggested by Bandelli and Rajagopal [18] : 

v = v(r,t)sin(9)e (p , (28) 

where (r, 9, ip) is a spherical coordinate system. Such an assumption can only be made if the 
unsteady part in the inertial term is neglected, or the inertial term is neglected as a whole. 

We consider the unsteady flow of the reacting mixture which is confined to the space between 
two cylinders, i.e., the mixture domain Q is defined as: ^ r ^ r Q , ^ 9 < 2n, and — oo < z < 
oo (see Figure 1). Here (r,9,z) represents a cylindrical-polar coordinate system that is placed 
such that r = rj is the boundary of the inner cylinder, which is fixed, and r = r is the inner wall 
of the outer cylinder (annulus), which is free to rotate. We shall employ a semi-inverse approach 
and seek a velocity field of the form: 

v = v e (r,t)e e + v z (r,t)e z , (29) 

where e# and e z denote the base vectors in the 9 and z coordinate directions, respectively. Note 
that this assumption for the velocity field automatically satisfies the isochoricity constraint (8). 
On substituting the similarity transformation (29) into the balance of linear momentum (23), 
expressed in a cylindrical-polar coordinate system, we obtain: 
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where 



ld_ 

r dr 



rfi c, 



dvg dv z \ dv 



dr ' dr J dr 



dir 



M I == *>(<=>< + 



dve_ _ve\ | ( dv z 
dr r J \ dr 



n(c) 



(30c) 



(30d) 



Note that in deriving (30), we have tacitly assumed that the pressure is independent of the 
9 coordinate. It would then be natural to assume, since we shall consider a fully-developed, 
oscillating flow, an axial pressure gradient of the form: 



^ = -[A + Bcos(f z t)], 



(31) 



where A and B are constants, and f z is the frequency of the oscillation. On integrating equation 
(31) we find the pressure: 



where h(r, t) satisfies: 



7r(r, z,t) — —[A+ Bcos(f z t)]z + h(r, t), 



dh = ^ 
dr r 01 



(32) 



(33) 



in view of the radial component of the balance of linear momentum (30a). If the circumferential 
component of velocity is known, then h(r, t) can be calculated via the formula: 



h(r,t)=g f riM ds+a ( t)) 

Jn s 



(34) 



where a(t) is an arbitrary function of time. The remainder of this paper is focused on determining 
Vg, v z , and c. To do so, we require both initial conditions and boundary conditions for vg, v z , 
and c; let us first consider the conditions on the velocity. 




Ug(t) 



Figure 1 . Illustration of the geometry for the problem under consideration along 
with a cylindrical-polar coordinate system. 
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5.2. Initial Conditions and Boundary Conditions on the Velocity Field. We shall as- 
sume that the fluid is initially at rest, and therefore the non-zero components of the velocity field 
must satisfy: 

v e (r, 0) = v z (r, 0) = VrG^rJ. (35) 

For fluids such as water flowing under "normal conditions," it is customary to enforce the "no- 
slip" condition! at the boundaries. However, if the mixture is sufficiently dense in either plasma 
or HA, then generalizations of Navier's slip condition may be more appropriate. For example, 
following the approach recommended by Rajagopal [19], one could have a boundary condition of 
the form: 

rw\T 

v rel ■ r + M([T T n • t\, c) U ' T = 0, (36) 

|T n • t\ 

where v re ; is the relative velocity between the boundary surface and the fluid, and r and n are 
the unit tangent and unit normal vectors, respectively, to the boundary surface. Note that M 
has been modified to account for the concentration of HA. In addition to the condition (36), the 
usual condition v ■ n = should also be satisfied if one assumes that the velocity of the mass 
diffusing through the articular cartilage and synovial membrane is negligible with respect to the 
velocity of the mass in the cavity. One example of the function M that has relevance for the 
problem under consideration is: 

M(\T T n.T\,c)= 1 —^\T T n.r\, (37) 

and one could require A(0) to be small so that when the fluid mixture primarily consists of 
plasma, the fluid nearly slips at the boundary. Also, if one requires A(l) = 1, then when the 
mixture is dense with HA, the no-slip boundary condition is enforced. With that being said, 
if the mixture is dense enough with HA, then the model represented by equation (20) may be 
ineffective in describing the response of the synovial fluid as viscoelastic effects could play a more 
dominant role. As a first order approximation we shall choose to enforce the no-slip boundary 
condition at the walls by requiring: 

ve{n,t) = and v s (r Q ,t) = r u g (t), t ^ 0, (38) 

where ug is the angular velocity of the outer cylinder. We shall assume that ug takes the form: 
LUg(t) = ujg[l — cos( fgt)], which could represent some cyclic routine such as exercise. Also, the 
cylinders are not permitted to translate, and therefore we must enforce the following boundary 
conditions for the axial component of the velocity: 

v z {r u t) =v z {r ,t) =0, t^0. (39) 

The solution to equations (30b) and (30c) ultimately depends upon the concentration, and this 
is the issue we take up next. 



§ While this condition is supposed to have had the backing of Stokes, he was far from unequivocal about its 
applicability to all flows of fluid such as water; he only advocated the same for sufficiently slow flows (see 
Rajagopal [191 for a detailed discussion of the same). 
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5.3. Assumptions for the Concentration. Since the mixture is confined to the annulur region 
between the two cylinders, and one of the cylinders is oscillating, we assume the concentration 
has the following structure: 

c=c(r,t), (40) 
and therefore the governing equation (26) reduces to: 

D c d ( dc\ dc 



dr\ T dr) 8f ^ 

5.4. Initial Conditions and Boundary Conditions for the Concentration. We shall re- 
quire the initial concentration be constant c» over the entire domain. The synovial membrane 
lines the joint cavity with the exception of the joint itself which is surrounded by articular carti- 
lage. If we take the inner cylinder to represent the joint, and the outer cylinder to be the cavity 
wall lined with synovial membrane, then the concentration must satisfy the following equations: 

dc 

— (n,t) = and c(r ,t) = Cr (t), t^O. (42) 

Condition (42) presumes that there is no diffusion of HA through the articular cartilage, which 
does in fact act as an HA reservoir. A more appropriate condition at the outer boundary would 
be to prescribe the concentration, until the mean value of HA within a region near the outer 
boundary reaches an optimum value at which time no more HA is secreted. This boundary 
condition would take the form: 

c{Toi^) 1 J C, C mean < C / .„ n 

— (r t) II c > c ' ^ > 

where 



Cmean(t) ■= =/ c(r,t)dr, (43b) 

r Q -r J- 

c is a constant, and c represents the optimal concentration level of HA within the layer r Q — f 
adjacent to the synovial membrane. Obviously the boundary condition (43) is a complicated 
condition to satisfy numerically as one must know, at a given instant of time, the entire con- 
centration profile. Thus, in addition to the iterations required due to the non-linearity of the 
equations, one must develop a iterative stencil to ensure that the condition (43) is satisfied, even 
if approximately. Thus, we shall first consider the more simple case wherein the conditions given 
by (42) must hold. 

5.5. Non-dimensionalization. We now proceed to non-dimensionalize the governing equations 
by introducing the mapping (r, z, t) >— > (f, z, t) and parameters: 

r-U Z Ir. V0 V z 7T 

t=f 6 t, V = ^^, W=^^, 7T 



r -ri V ' r o 0Je r o 0Je Q f r uo e Lf e ' 

c = c{r,t) = c(r,t), /u = — , u e = — , Pf = ~r, Vg = , (44) 

u e J e r ~ri 

2jr 2 u 2 e 
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where L is a characteristic length. First, we record the non-dimensional solution for the pressure 
for completeness: 



7t(r, z, t) = -[p A + p B cos(p f i)]z + 



r -\2 



v 2 (s,t) 
o *+Pg 



ds + a(t), 



(45) 



where a(t) :- 



q r'weLfg i anc ^ ^ e P ressure gradient coefficients are given by the formulae: 



.4 



PA 



and pb 



B 



Q f r u e fe ' QjroUJefe 

The non-dimensional forms of the governing equations (30b), (30c), and (41) are: 



1 







Re (f + pg) 2 dr 



K r +Pg) V c , 



dv dw\ ( dv 



dr' dr J \dr f + p g 



dv 



and 







Re (f + p g ) df 



(r + Pg)P> c 



dv dw\ dw 



dr' dr ) dr 



1 1 d 

Pe f + p g dr 



[r + p g ] 



dc 
dr 



+ PA+PBCOs(p f t) 

dc 



dw 



respectively, where 



dv dw 



dv v 



2 ( dw~ 1 
\ df 



dr f + pg 

and the Reynold's and Peclet numbers are defined as follows: 

Re:= gMr °~ ri)2 and Pe := " ^ 



n(c) 



/'() 



D, 



(46) 

(47a) 

(47b) 
(47c) 

(47d) 
(47e) 



respectively. Below we summarize both the initial conditions and boundary conditions, in non- 
dimensional form, for v and w: 



v(f,0) = w(r,0) = 0, V f E [0, 1], 

v(0, i) = and v(l, t) = 1 - cos(t), t ^ 0, 
w(0,i) = w(l,i) = 0, t^0, 



and for the concentration c: 



9c 



c(f,0) = 0.1, V f G [0,1], 



— (0,t) = 0, t ^ and c(l,t) 



0.1 + 0. It, fe [0,2] 
0.3, £ > 2 



(48a) 
(48b) 
(48c) 

(49a) 
(49b) 



Furthermore, we shall require pa = ~Pb, and ps = 1, and these parameter values agree with the 
initial condition (48a). The ramp function given in condition (49b) is the same used by Hron et 
al. [S] and this, to some extent, represents physiological data. 
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6. Discussion 

Equations (47a)-(47c) constitute a non-linear system of partial differential equations that must 
be solved subject to the previously discussed initial conditions and boundary conditions, for both 
the velocity and the concentration. We use the subroutine pdepe within the software package 
MATLAB. This package utilizes a Petrov-Galerkin finite element method to generate a system 
of ordinary differential equations [20] . The ordinary differential equations are then solved using 
standard methods. 

The set of coupled partial differential equations are solved for three shear-thinning, chemically- 
thickening models (Model 1, Model 2a, and Model 2b), and the results are compared with those 
for the Newtonian model. Note that since the diffusivity was assumed to be constant, the 
concentration profiles are the same for all the models (see Figure (13d)) ) . i.e., while the components 
of velocity depend upon the concentration, the converse is not true. This is because the equation 
governing concentration is not influenced by the velocity field as the diffusivity is a constant, 
and a special form is assumed for the concentration and the velocity field. If one allows for the 
diffusivity to depend on the norm of the first Rivlin-Ericksen tensor, |Ai|, then the equations 
would be further coupled, and the concentration profiles would vary depending on the model, 
at any given instant of time. The concentration profiles will also further change if an external 
pressure gradient is applied in the z direction, when one assumes the concentration to depend 
on the norm of the first Rivlin-Ericksen tensor. Furthermore, note that one cycle is complete 
in (non-dimensional) time of 2n, regardless of whether or not there is an axial pressure gradient 
present. 

We shall first consider the case wherein there is no pressure gradient applied in the z direction. 
The non-dimensional velocity profiles in the 6 direction are depicted in Figures (|2J), (13"a|) . and 
f!3bj) . The non-dimensional velocity in the z direction is zero for all cycles for each of the models 
studied because there is no pressure gradient in the z direction, and the axial velocity is zero at 
the boundaries. The non-dimensional apparent viscosity variations are compared in Figures (j4l 
and dSJ). Obviously the plots for the chemically-reacting fluid illustrate a great departure from 
those associated with the Newtonian model, and this is to be expected in view of the change in 
the concentration over the domain (Figure (13d[) ) . Now, as the number of cycles increases, more 
reactant enters the domain from the outer boundary, and this reactant stays in the domain as 
the inner cylinder is non-porous. This, of course, increases the concentration of the reactant in 
the region between the cylinders (see Figure (13d)) ) . which in turn increases the apparent viscosity 
as observed in Figures (jlj) and (jSJ). Hence, as the cycles continue, the region near the outer 
boundary becomes more and more "viscous." This outer layer of highly viscous fluid squeezes 
the less viscous fluid particles, radially inward, and causes them to accelerate in the 8 direction. 
Such behavior is illustrated in Figures (j2J), fl3a|) . and f!3bj) . wherein the velocity in the 9 direction 
increases in the central region with number of cycles. Such phenomena was also observed by 
Bridges and Rajagopal [8 J while studying shear-thinning/chemical-thickening fluids. However, 
for their model, the zero shear-rate viscosity had a quadratic dependence on the concentration, 
and the shear- index was assumed constant. Furthermore, it was assumed that the diffusivity 
depends on the first Rivlin-Ericksen tensor. In addition, the geometry and boundary conditions 
were slightly different than those for the problem under consideration. 

Next, we shall study the case when there is an external pressure gradient in the z direction. 
The non-dimensional velocity profiles in the 6 direction are shown in Figures ([6]) and (J7|), and the 
non-dimensional velocity profiles in the z direction are compared in Figures jS]), (19a)) . and (19bj) 
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for the various models under consideration. The non-dimensional apparent viscosities in this case 
are compared in Figures (fit)]) and (III]) . One can also infer from Figures (|HJ), (|9aj) . (j9b]l . and ( )9d]) . 
that as the number of cycles increase, the component of velocity in the z direction decreases in the 
central region. Now, w is zero at the boundaries, and is driven solely by the pressure gradient. 
As the number of cycles increases, there is a significant increase in the apparent viscosity in 
the central region. This increase in viscosity greatly mitigates the driving force of the pressure 
gradient, and thus decreases the velocity w. Furthermore, in comparison with the results from 
the case when there is no pressure gradient, it can be seen that at fewer cycles, there is significant 
difference in the velocity profiles in the 9 direction and the apparent viscosity. For instance, one 
can see this difference by comparing Figures (I2dj) and (I6dj) at 3.5 cycles. However, for more 
cycles (12.5, 34.5), the plots look similar. The reason being that the velocity in the z direction, 
and its gradients are steep, at lower cycles. This phenomena affects the apparent viscosity, and 
hence in turn the velocity in the 9 direction. At higher cycles, as discussed above the velocity in 
the z direction is small, and its gradients are small as well; thus, there is no significant effect on 
the apparent viscosity and the velocity in the 9 direction. 

When one compares the three shear-thinning, chemically-thickening models, the velocity and 
the apparent viscosity profiles for Model 2a and Model 2b are similar, and those of Model 1 
are quite distinct from Models 2a, 2b; this was something Hron et al. [6] also observed. It is 
also interesting to note that the nature of the shapes for the velocity plots for shear-thinning, 
chemically-thickening models, and the Newtonian model is completely different (see Figures fl2]), 
( )3a| . and (13bp . for instance) - the plots for Models 1, 2a, 2b are concave and the plots for the 
Newtonian model are convex. It is quite possible that this is due to the curvature in the flow 
domain since the results for a rectangular domain as seen in [B] do not show such a behavior. In 
the study by Hron et al. [S] , they find that the structure of the solutions for all the solutions for 
the three models are similar. Thus, studying the fluid in a different geometry seems to provide 
some insight in that the structure of the solution for the non-Newtonian fluids are markedly 
different from those of the Newtonian fluid. We were hoping that the curvature effects would 
help in providing us with new information concerning the response of these non- Newtonian fluid 
models and our guess was correct. These non-Newtonian fluids are quite different from the 
Newtonian fluid as evidenced by the profiles in the non-Newtonian case being concave and the 
Newtonian case being convex. This implies that the nature of the shear stresses in the non- 
Newtonian case and the Newtonian case are totally different. It would be useful to find an 
initial boundary value problem where the three different non-Newtonian models lead to different 
solutions with the hope that one agrees better with the corresponding experimental results than 
the two others, thereby allowing us to choose an appropriate model. 

This study is merely a preliminary study to assess a fluid model that can describe the response 
of fluids such as the synovial fluid. This study will be followed by an analysis of such a fluid in a 
realistic geometry, the resolution of which will involve numerical resolution based on appealing 
to methods such as the finite element method. 
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Figure 2. Non-dimensional velocity profiles in the 9 direction. Comparison for 
different models at different cycles are shown in (a), (b), (c). In (d), the non- 
dimensional velocity profiles in the 9 direction are compared for different cycles 
for Model 1. The parameters chosen are = 1, r Q = 1.2, pa = ~Pb — 0, p/ = 1, 
Uq — 1, Re — 10, Pe = 1000. The pressure gradient in the z direction is chosen to 
be zero here. 
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Figure 3. The non-dimensional velocity profiles in the 9 direction are compared 
for different cycles. Results for Model 2a and Model 2b are shown in (a) and (b) 
respectively. Non-dimensional centerline velocity (at f = 0.5) in the 9 direction as 
a function of number of cycles for different models is shown in (c). Concentration 
profiles at different cycles are shown in (d). Since, the diffusivity is chosen to 
be a constant, the concentration profiles are the same for all the models. The 
parameters chosen are = 1, r a = 1.2, pa = ~Pb — 0, pf = 1, oJq = I, Re = 10, 
Pe = 1000. The pressure gradient in the z direction is chosen to be zero here. 
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(c) 34.5 cycles 

Figure 4. Non-dimensional apparent viscosity profile comparison for different 
models at different cycles are shown in (a), (b), (c). The parameters chosen are 
n = 1, r = 1.2, p A = -p B = 0, p f = 1, u e = 1, Re = 10, Pe = 1000. The 
pressure gradient in the z direction is chosen to be zero here. 
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Figure 5. The non-dimensional apparent viscosity profiles are compared for dif- 
ferent cycles for Model 1, Model 2a and Model 2b. Results are shown in (a), (b) 
and (c) respectively. Non-dimensional centerline apparent viscosity (at f = 0.5) as 
a function of number of cycles for different models is shown in (d). The parameters 
chosen are = 1, r D = 1.2, p^ = ~Pb — 0, Pf — 1, uj$ = 1, Re = 10, Pe = 1000. 
The pressure gradient in the z direction is chosen to be zero here. 
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Figure 6. Non-dimensional velocity profiles in the 9 direction. Comparison for 
different models at different cycles are shown in (a), (b), (c). In (d), the non- 
dimensional velocity profiles in the 9 direction are compared for different cycles 
for Model 1. The parameters chosen are = 1, r Q = 1.2, pa = ~Pb — 1> Pf — 1, 
us = 1, Re = 10, Pe = 1000. The concentation profiles for this set of results is 
same as (l3d|) because the diffusivity is chosen to be constant. 
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(a) Comparison for Model 2a 
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Figure 7. The non-dimensional velocity profiles in the 6 direction are compared 
for different cycles for Model 2a and Model 2b are shown in (a) and (b) respectively. 
The parameters chosen are = 1, r Q — 1.2, pa = ~Pb — 1> Vf — 1> — 1> 
Re = 10, Pe = 1000. 
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Figure 8. Non-dimensional velocity profiles in the z direction. Comparison for 
different models at different cycles are shown in (a), (b), (c). In (d), the non- 
dimensional velocity profiles in the z direction are compared for different cycles 
for Model 1. The parameters chosen are = 1, r Q = 1.2, pa = ~Vb — 1, Pf = 1, 
uj e = 1, Re = 10, Pe = 1000. 
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Figure 9. The non-dimensional velocity profiles in the z direction are compared 
for different cycles for Model 2a and Model 2b. Results are shown in (a) and (b) 
respectively. Non-dimensional centerline velocity (at f = 0.5) in the 9 direction as 
a function of number of cycles for different models is shown in (c). Non-dimensional 
centerline velocity in the z direction as a function of number of cycles for different 
models is shown in (d). The parameters chosen are = 1, r a = 1.2, p^ = ~Pb = 1, 
p f = l,u g = 1, Re = 10, Pe = 1000. 
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(c) 34.5 cycles 

Figure 10. Non-dimensional apparent viscosity profile comparison for different 
models at different cycles are shown in (a), (b), (c). The parameters chosen are 
n = l,r Q = 1.2, pa = -Pb = l,Pf = 1, u e = 1, Re = 10, Pe = 1000. 



ON MODELING THE RESPONSE OF SYNOVIAL FLUID 



25 



7 1 1 1 1 6 




— 1 1 1 1 1 1 1 — — 1 — — 1 — — 1 — — 1 — 

0.2 0.4 0.6 0.8 1 0.2 0.4 0.6 0.8 1 



r r 
(a) Comparison for Model 1 (b) Comparison for Model 2a 




f No. of cycles 

(c) Comparison for Model 2b (d) Centerline fi for different models 



Figure 11. The non-dimensional apparent viscosity profiles are compared for 
different cycles for Model 1, Model 2a and Model 2b. Results are shown in (a), (b) 
and (c) respectively. Non-dimensional centerline apparent viscosity (at f = 0.5) as 
a function of number of cycles for different models is shown in (d). The parameters 
chosen are r, = 1, r a = 1.2, p A = —p B = 1, pf = 1, u e — 1, Re — 10, Pe = 1000. 
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